
#include <iostream>
#include "matrix.h"
#include "matrix_operations.h"
#include "matrix_compound_operations.h"
#include "eigenvalue.h"

// generelle Hinweise zum LA-Paket:
// 1. es gibt folgende Matrizenklassen:
//        Matrix (allgemeine Matrix, abstrakt, nur lesbar)
//        ReMatrix (allgemeine Matrix, les- und schreibbar, von Matrix abgeleitet)
//        ColumnVector (Spaltenvektor, les- und schreibbar, von ReMatrix abgeleitet)
//        RowVector (Zeilevektor, les- und schreibbar, von ReMatrix abgeleitet)
//        SymmetricMatrix (symmetrische Matrix, abstrakt, nur lesbar)
//        SyMatrix (symmetrische Matrix, les- und schreibbar, von Matrix und SymmetricMatrix abgeleitet)
//        UpperTriangularMatrix (obere Dreiecksmatrix, abstrakt, nur lesbar)
//        UTMatrix (obere Dreiecksmatrix, les- und schreibbar, von Matrix und UpperTriangularMatrix abgeleitet)
//        LowerTriangularMatrix (untere Dreiecksmatrix, abstrakt, nur lesbar)
//        LTMatrix (untere Dreiecksmatrix, les- und schreibbar, von Matrix und LowerTriangularMatrix abgeleitet)
//        DiagonalMatrix (Diagonalmatrix, abstrakt, nur lesbar)
//        DiMatrix (Diagonalmatrix, les- und schreibbar, von Matrix, SymmetricMatrix, DiMatrix, 
//                  UpperTriagularMatrix, LowerTriangularMatrix abgeleitet)
//
// 2. Zugriff auf die Matrizen:
//    folgende Methoden stehen zur Verfuegung:
//      - number_rows(), number_columns() liefern Zeilen/Spaltenzahl
//      - clear()  setzt alle Matrixeintraege auf 0 zurueck
//      - entry(i,j) liefert den Eintrag in der i-ten Zeile und j-ten Spalte (Numerierung ab 1 !!!!),
//                   bei ColumnVector und RowVector: entry(i)
//      - r_entry(i,j) liefert eine Referenz auf den Eintrag in der i-ten Zeile und j-ten Spalte (falls moeglich),
//                     bei ColumnVector und RowVector: r_entry(i)
//                     Zugriff auf die "0"-Eintraege in UTMatrix, LTMatrix, DiMatrix mit r_entry fuehrt zu Problemen !!!!
//      - at(i,j), r_at(i,j) tun das selbe wie entry(i,j), r_entry(i,j), pruefen allerdings, ob die Indeces im
//                           erlaubten Bereich sind und werfen ggf. eine Ausnahme
//
// 3. in den Headern finden sich:
//    "matrix.h": die Deklarationen der Matrixklassen, keine Operationen
//    "matrix_operations.h": einfache Operationen wie Kopieren, Addieren, Subtrahieren, Multiplizieren, 
//                           Skalarmultiplikation, Skalarprodukt, Inverse, Cholsky-Zerlegung, Transponieren, 
//                           Determinante, Spur
//    "matrix_compound_operations.h": zusammengesetzt Operationen, z.B. Multiplikation mit Transponierter,
//                                    mehrfache Addition u.v.m.
//    "matrix_norms.h": einige Vektornormen und Abstandsmase (L1,L2,Maximum,Minkowski)
//    "eigenvalue.h": Eigenvektor/Eigenwert-Berechnung fuer positiv definite Matrizen
//
// 4. Die Matrixoperationen besitzen ueblicherweise den Aufbau:
//    void <Funktionsname> ( <Zielmatrix>, <Argument1>, <Argument2>, ...)   wenn das Ergebnis eine Matrix ist
//    double <Funktionsname> ( <Argument1>, <Argument2>, ...)   wenn das Ergebnis ein Skalar ist
//    in "matrix_compond_operations.h" sind einige Funktionen deklariert, bei denen das erste Argument 
//    sowohl Eingabe- als auch Zielwert der Funktion ist, z.B. add_to, add_to_multiply_with_transpose, u.v.m
//
// 5. Fehlerfreiheit kann natuerlich nicht garantiert werden.



using namespace Algebra;   // Namespace bekannt machen
using namespace std;


int main () {
  try{

    // einfach mal einen 4-dimensionalen Spaltenvektor
    ColumnVector vec1 (4);
    // die Eintraege setzen: Indizierung von 1..4
    for (unsigned int i=1; i<=vec1.number_rows(); i++)
      vec1.r_entry(i) = 0.3*static_cast<double>(i);  // r_entry liefert eine Referenz auf einen Vektor/Matrixeintrag
    
    // zur Abwechslung mal eine 3x4 Matrix (3 Zeilen, 4 Spalten):
    ReMatrix mx1 (3,4);   // "Re" steht fuer rectangular (rechteckig)
    for (unsigned int i=1; i<=mx1.number_rows(); i++)
      for (unsigned int j=1; j<=mx1.number_columns(); j++)
	mx1.r_entry(i,j)=(i*j)%7;

    // Multiplizieren:
    ColumnVector vec2 (3);
    multiply (vec2, mx1, vec1);   // standardmaessige Argumentfolge bei einfachen Operationen: Ergebnismatrix, Argument1, Argument2, (...)
    
    // eine symmetrische Matrix konstruieren durch Multiplikation von vec2 mit vec2/transponiert:
    SyMatrix mx2 (3);   // 3x3 symmetrische Matrix
    multiply_with_transpose (mx2, vec2);   // in "matrix_compound_operations.h" sind verschiedene Kombinationen einfacher Operationen definiert, z. B. Multiplikation mit der Transponierten, siehe den Header

    // Die Einheitsmatrix erzeugen:
    DiMatrix id (3);
    identity (id);

    // id zu mx2 addieren
    add_to (mx2, id);

    double d = determinant (mx2);   // Determinante berechnen
    double s = trace (mx2);         // Spur berechnen

    // Eigenwerte und Eigenvektoren der Matrix mx2 bestimmen (bisher nur fuer symmetrische Matrizen implementiert)
    ColumnVector vec3 (3);
    ReMatrix mx3 (3,3);
    eigenvalue (vec3, mx3, mx2);   // vec3 enthaelt die Eigenwerte von mx2, mx3 spaltenweise die zugehoerigen Eigenvektoren

    // Ausgeben der Ergebnisse:
    cout << "Determinante: " << d << '\n';
    cout << "Spur: " << s << '\n';
    ColumnVector vec4 (mx3.number_rows());
    for (unsigned int i=1; i<=vec3.number_rows(); i++) {
      for (unsigned int j=1; j<=mx3.number_rows(); j++)
	vec4.r_entry(j)=mx3.entry(i,j);    // kopieren der i-ten Spalte von mx3 nach vec4
      cout << i << ". Eigenwert: " << vec3.entry(i) << " mit Eigenvektor " << vec4 << '\n';
    }
  }catch(exception& e){
    // falls irgendwas schief gegangen ist, zum Beipiel wenn eine Matrix falscher Groesse uebergeben wurde, fange die geworfene Ausnahme
    cerr << e.what() << endl;
    return -1;
  }
  return 0;
}
